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Abstract: We present a new tree-level matrix element generator, based on the colour 
dressed Berends-Giele recursive relations. We discuss two new algorithms for 
J> ' phase space integration, dedicated to be used with large multiplicities and 

O . colour sampling. 
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1 Introduction 

43 : 

In recent years considerable progress has been made in the calculation of full matrix elements (ME) for higher 
order perturbative corrections to Standard Model (SM) processes, QCD and QCD associated processes in 
particular. Automatic computation of NLO virtual corrections to arbitrary processes finally seems within 
reach due to newly emerging numerical techniques [TJ2] . On-shell recursive methods proved to yield compact 
expressions for multi-leg tree-level amplitudes with massless [3J and massive [I] external particles and are 
CNl ■ now widely used. The CSW vertex rules [5] as off-shell techniques are employed in many analytical and 

numerical approaches 

Apart from major developments in the computation of loop amplitudes, many attempts have been made to 
tackle the task of numerically evaluating tree-level amplitudes with large numbers of external legs. They 
led to the construction of several programs, capable of evaluating general tree-level processes [c^l^ [TOI[TT] . 
In this context it turned out, that with increasing number of particles involved in the scattering one of the 
the most efficient methods to compute colour-ordered amplitudes is the Berends-Giele recursion [T2 j [T3 l [T4 ] . 
Correspondingly the fastest methods available for the computation of full scattering amplitudes are the 
colour dressed Berends-Giele relations [TS] , which are essentially equivalent to the Dyson-Schwinger methods 
employed in Refs. [To] , with the ALPHA algorithm of Ref. [T7] being comparable in efficiency. In Refs. [H] 
and [E] it was pointed out that a vertex decomposition of four-gluon vertices in QCD is clearly advantageous 
if the speed of numerical implementations is concerned. These findings raise the question, whether it is 
possible to construct a full set of SM Feynman rules with no four vertices present in the theory, such that 
recursive techniques analogous to the colour dressed Berends-Giele relations can be employed in numerical 
programs. In Sec. [5] we demonstrate that this is feasible. We discuss the numerical implementation of the 
results in the new ME generator COMIX in Sec. [3] and present code-related aspects, such as a multi-threading 
concept. 

A very important part of computing cross sections for tree-level processes is, to find an efficient algorithm 
for phase space generation. If colours are sampled over, similar problems arise for colour space. An effective 
general technique for phase space generation has been presented in Ref. [TH] . We observe in Sec. 14.11 that 
it is possible to formulate the rules presented ibidem in a truly recursive fashion, i.e. on the same footing 
as the matrix element computation. This implies in particular, that point by point the same calculational 
effort is spent for computing matrix clement and phase space weight. We introduce effective colour sampling 
techniques in Sec. 14.21 Having these techniques at hand, we elaborate on how to eventually couple colour and 
phase space integration and propose a new type of integrator based on the HAAG generator [19j in Sec. [ 
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* See |http: / / comix. freacafe.de| for downloads and a detailed manual. 
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We present a comprehensive comparison of results generated with COMIX to those generated with the two 
other multi-leg tree-level matrix element generators AMEGIC++ [H] and ALPGEN [TT] in Sec. [5] Section \§\ 
contains our conclusions. 



2 Recursive relations for tree-level amplitudes in the Standard Model 

It has been pointed out, for example in Refs. [IH[T4l[T5], that the calculation of multi-parton amplitudes 
is substantially simplified when employing Berends-Giele type recursive relations. One main reason for 
the simplification is that these relations allow to reuse basic building blocks of an amplitude, which are 
the m-particle internal off-shell currents. Another reason is that they can be easily rewritten to include 
three-particle vertices only. In the following we will briefly illuminate, why this is a major advantage. 



2.1 The cost of computing a tree amplitude 

As an example, we try to estimate the total computational cost for tree amplitudes, given a certain type of 
vertices in the underlying theory. We assume that only one particle type exists and the internal n-particle 
currents obey a recursion, which is of the functional form 

n 

J n (it) = P n (n) ^2 X! (ki,---,tn) Jii (tti).-. Ji N (ttjv) ■ (1) 

iV=lPjv(7r) 

Here J m denote unordered m-particle currents, while Vn are N + 1-point vertices and P n is a propagator 
term. The two sums run over all possible vertex types Vn and all (unordered) partitions Vn (t) of the set 
of particles tt into N (unordered) subsets, respectively [IS]. The full n + 1-particle scattering amplitude 
can be constructed by putting an arbitrary n-particle internal off-shell current on-shell and contracting the 
remaining quantity with the corresponding external one-particle current. 

A n+1 {tt) = Ji (t) p ^ * J n (tt \ i) . (2) 

We now deal only with vertices of N + 1 external legs and we consider their contribution to the computation 
of an n-particlc off-shell current. The number of vertices to evaluate per m-particle subcurrent is the Stirling 
number of the second kind S (to, N), corresponding to the number of partitions of a set tt of m integers into 
N subsets. The total number V(n, N) of N + 1-particle vertices to be calculated thus becomes 



m—N 



Since the Stirling numbers S{m, N) are zero for m < N, we can extend the sum down to zero, leading to 



'N 



rn=0 v 7 
1 

~ (N + l) 



N 



(N - i) m 

(4) 

^EKf^ 1 ) {N+l-i) n+1 =S(n + l,N + l) . 



i=0 



The question is, whether we can obtain a milder growth in computational complexity, if all N + 1-particle 
vertices occuring in Eq. (fT]) are decomposed in terms of two or more vertices with fewer number of external 
legs. When doing so, we must introduce additional pseudoparticles reflecting the structure of the decomposed 
vertex. Hence we have to consider the contribution arising from the presence of these pseudoparticles, too. 
The problem can be simplified by assuming that there is only one additional pseudoparticle, which obeys a 
completely independent recursion relation. Then the full contribution of an N + 1-particle vertex, now being 
decomposed into a M + 1- and a N — M + 1-particle vertex becomes 

S (n + 1, N + 1) -> S (n + 1, M + 1) + S (n + 1, N - M + 1) , (5) 
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which can be either bigger or smaller than S (n + 1, N + 1), depending on n, N and M. With increas- 
ing n, however the right hand side is always smaller such that the vertex decomposition becomes clearly 
advantageous. Similar arguments hold when introducing more than one pseudoparticle. 

From this simple but general consideration we see that the aim of any recursive formulation of interaction 
models should be, to reduce the number of external lines at interaction vertices to the lowest possible. In 
this section we will show that within the Standard Model it is possible to reduce iV max to two, which is the 
lowest possible number in general. For QCD interactions we employ the results of Ref. |15j . where this task 
has already been performed and the original Berends-Giele recursive relations have been reformulated to 
incorporate colour. 



2.2 General form of the recursive relations 

In the following we will denote by J a (it) an unordered SM current of type a, which receives contributions 
from all Feynman graphs having as external particles the on-shell SM particles in the set tt and one internal 
particle, described by this current. The index a is a multi-index, carrying information on all quantum 
numbers and eventually on the pseudoparticle character of the particle. Special currents are given by the 
external particle currents. They correspond to external scalars, spinors and polarisation vectors, see Sec. [3] 
For them there is only one multi- index a = a% associated with the external particle i, whereas in the general 
case multiple multi-indices may lead to non- vanishing internal currents. This corresponds to multiple particle 
types being possible as intermediate states. Assuming that only three-point vertices exist, any internal SM 
particle and pseudoparticle off-shell current can be written as 

Ja (7T) = P a (7T) J2 S fa' 71 *) K? 1,0,3 fa' 71 *) J ^ fa) J ^ fa) • ( 6 ) 

V™ 1,0,2 PaO) 

Here P a (tt) denotes a propagator term depending on the particle type a and the set tt. The term V^ 1 '™ 2 (tti, tt-i) 
is a vertex depending on the particle types a, ot\ and a.2 and the decomposition of the set tt into disjoint 
subsets 7Ti and tt2- The quantity S (tti, 7^) is the symmetry factor associated with the decomposition of tt 
into 7Ti and 7T2 and will be discussed in Sec. 12.51 Superscripts in this context refer to incoming particles, 
subscripts to outgoing particles. The sums run over all vertices in the reformulated Standard Model and 
all unordered partitions V2 of the set 7r into two disjoint subsets, respectively. A full unordered n-particle 
scattering amplitude is then given by 

A (tt) = J an (n) p _ ^ J an (tt \ n) , (7) 

where a denotes a set of reversed particle properties, i.e. opposite helicity, colour, momentum and particle 
type. It has been proved in Ref. [15] that the above form is correct for pure gluonic scattering amplitudes 
once the four gluon vertex is suitably decomposed into two vertices involving an internal antisymmetric 
tensor pseudoparticle. We briefly recall this proof before continuing with the decomposition of four particle 
vertices in electroweak interactions. Once this decomposition is achieved, no further complications arise and 
Eq. © can be employed to compute arbitrary scattering amplitudes in the Standard Model. 



2.3 Colour dressed Berends-Giele recursive relations in QCD 

Any perturbative QCD scattering amplitude A can be written as a sum of terms, which factorise into two 
components, one only depending on the gauge structure and one only depending on the kinematics. Such 
a decomposition is called colour decomposition. Considering for example tree-level n-gluon amplitudes, 
several colour decompositions exist. A very intuitive one based on the fundamental representation of the 
gauge group is given by [20] 

A(l,...,n)= Tr (T ai T a °* ...T Q -) A (1,(72,.-., <r n ) ■ (8) 

<?es„_i 

Here a runs over all permutations of the n— 1 indices 2 ... n. The functions A depend on the Lorentz- 
structure of the process only and are called colour-ordered amplitudes. A more suitable colour decomposition 
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for n-gluon amplitudes has been introduced in Refs. [2"T] . It employs the adjoint representation matrices 
(F a ) bc of SU(3) and reads 

A(l,...,n)= J2 ( pa ° 2 ■■■ Fa ° n - 1 )a ian A(l,a 2 ,...,a n ^,n) . (9) 

Note that in this case the sum runs over the permutations of the n — 2 indices 2 ... n — 1 only, whereas the 
first and the last index remain fixed. Another colour decomposition, suited especially for Monte Carlo event 
generation is the colour flow decomposition [22] . In this prescription the SU (3) gluon field is treated as a 
3x3 matrix (A^) 1 ^ rather than a one index field A". The corresponding decomposition reads 

A(l,...,n)= S ll ^5 1 '^ ...S^nh A(l,a 2 ,...,a n ) . (10) 

<?es„_i 

The remaining task is now, to compute the colour-ordered amplitudes. In Ref. [T^] Berends and Giele 
proposed a method to do so in a recursive fashion. The basic idea is that, according to the Feynman rules 
of QCD, an internal n-gluon current is defined by all contributing Feynman graphs with n external on-shell 
gluons and one off-shell gluon. 

fn-l 
y; vr x ( Phk ,p k+1 . n ) j k a, . . . , *o Jx (t + 1, . . . , n ) 
k=i 

( n ) 

7i — 2 n—1 

J2 E v: pKX j p (i,...,j)j K (j + i,...,k)j x (k+i,. 

3 = 1 k=j+l 

Here pi denote the momenta of the gluons, pi.j = pi + ... + Pj and V£ kX (pi,k,Pk+i,n) and y^ pliX are the 
colour-ordered three and four-gluon vertices defined according to Ref. [13] , 

(jp,q)=i^=( 9 KX ( P -q) U + 9 Xu (2q + pf - g VK (2p + q) X ) , 

f ' (12) 

yUpKX = • ^ ( 2g „ KgP X _ g U PgK X _ gV X gP , ) 

The full colour-ordered n-gluon amplitude A (1, . . . , n) is then obtained by putting the n— 1-particle off-shell 
current J n -i (1,- •• 1) on-shell and contracting it with the external polarisation J M (n). Employing the 
tensor-gluon vertex 

^ KA = ^(ff M V A -^V K ) , (13) 

and the tensor "propagator" 

-W£ = -i(g&-g*gZ) , (14) 
the recursion can be reformulated to give 

n-l ( 

-T-g^v 

k=1 K \ (15) 



J M (1, 2, . . . , n) = \ V ™ X (Pi,*»P*+l.n) J K (l,...,k)J X (k + l,. 



+ V™ al3 J K {l,...,k)J af3 (k + l,...,n) + V XvaP J a p (1, . . . ,k) Jx {k + 1, . . . ,n) | 

and 

J af} (1, 2, . . . , n) = -iD«i E V^ kX J k (1, . . . , k) Jx (k + 1, • . . , n) , (16) 



k=l 



for the gluon and tensor pseudoparticle currents, respectively. Since no external tensor currents exist, all 
tensor currents with one particle index only are defined as zero. The advantage of the above formulation in- 
cluding a tensor current, as discussed in Sec. 12. H is the elimination of the four-gluon vertex. Correspondingly 
we introduce a "pseudogluon" , which, from here on, we denote by 174. 
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Following Ref. [TS], one can introduce colour dressed gluon and tensor pseudoparticle currents J^jj and 
Jap I J, defined by 

J»ij(l,---,n) = ^2 6 IJ<ii di CTi j iT2 ...d ian j J fl ((j 1 ,...,a n ) , 

(17) 

Denoting by ir the set (1, . . . , n) of n particles, the following recursive relations for these currents are obtained: 

, vKKL,apMN n / w / \ I (18) 

0-P 2 (t) 

J a pijfa = D 2l"^ E V^5 K HG XMN ^KKlfa) J\MN fa) ■ 

Here we have defined the colour dressed gluon and tensor pseudoparticle vertices 

K,ff<5 (7Ti,7r 2 ) = dgtf %T/ 3l/ (7ri,7T 2 ) +%d V 3l/ (7r 2 ,7Ti) , (19) 

and 

<S " M * - # + 6%6 ML SE V T % . (20) 

The second sum runs over the set of ordered partitions of the set ir into two disjoint subsets, OV2 fa- 
A complete proof of these relations can be found in Ref. [15] . The above procedure of colour dressing 
can easily be generalised to QCD processes including quarks. Since no further elementary QCD four-point 
interactions exists, no further vertex decomposition has to be performed and therefore no new current types 
are introduced. For amplitudes including quarks care must be taken of using the proper colour space gluon 
propagator when coupling to qqg vertices, i.e. 

P.j? « Sfsf-^SjjS^, (21) 
as described in Ref. [25] , 

2.4 Decomposition of electroweak four-particle vertices 

The above procedure can be generalised to describe all Standard Model interactions, once a suitable replace- 
ment of the corresponding four particle vertices has been found. 

We start by proposing a decomposition of four particle vertices with VK-bosons onljQ 

, ; W~k, W + v. W~\ , ; W~k,Z47<5 p <*P ^,W + v.W~\ , ^jW~ X, Z^S p a/3 , ; W+k,W~k /nni 
y W-fj, ^ V W~ il ' ' V Z ia ,p + V W~n ' r Z Al 8 ' V Z A a.p • V ZZ J 

Here denotes a new antisymmetric tensor pseudoparticle introduced for the vertex decomposition. Its 
interaction vertex reads 

v w-n ~2 ^ M M ' ' z i a P ~2 \9a9p 9 a 9p) ■ l ZJ 7 

To obtain correct signs of four-particle vertices, we define the tensor pseudoparticle "propagators" as 

P a $ = **D$ where if °£ ^ , (24) 



1 Note that this decomposition of vertices is not unique and other choices may exist. 



■5 



94 - %4 

fu/, IJ k\,HG = -iSfStjDft /'" ,,A = i D$ 



K A 



= i I'" ''A = iD^ 

Tab. 1 Standard Model propagators for auxiliary particles introduced in the vertex decompo- 
sition. Note that the l/Afc-term arising in Eq. (|2ip is obsolete for the pseudogluon 
propagator because the pseudogluon does not couple to quarks. 



and where is given by Eq. (fT4|) . Note that the Z 4 pseudoparticle is not self-conjugate. This definition 
prevents double counting four-particle vertices involving the W boson and constructing fake WWWW 
vertices with all W's having the same charge. The four-particle vertices involving W bosons, photons and 
Z-bosons are decomposed as follows 

^ AK. VV V . Zj A y ™, (» 4 70 „ up ->jW V,£* I \y 1" p ap j VV U,JiK (OK, 

V W-u. V W-u. ■ "-or- ~x ' V \r,r-r,n V W-u ' ~vt/- ~x ' V \iir-„a ' V Z0 J 

We introduced a new tensor pseudoparticle, , whose interaction vertices are defined as 



)Ak, W v, AX 


Ak,WZ 


76 




a/3 


.v w ' 


u,A\ 


V A\,W 4 ~ 


' ~y8 




a/3 


,v w ' 


v, Ak 






■yS 


a/3 




76 


a/3 


,Ak, W~v, Z\ 
W~n 


V W- i x 


7(5 




a/3 
7<5 


,v w 


ZA , 

a/3 


v z\,w; 


7(5 




a/3 
76 


■ V w ~ 


~ v , Ak 

a/3 


, Zk, W~u, ZX 


v Zk,WZ 


7(5 




a/3 


,v w 


~ v. ZX j 


' v w-n 


7<5 


■v 


a/3 


v w 


~ v, Zk, 


W-/J. 




7(5 


a/3 




7(5 


a/3 



.,Ak.W 4 7<5 * . - / ~ K § g .,w v,Am 1 ■ a (uk k u\ 

V W -^ = 2^ sm0 w [g^g -g^g 1 ), v Wi - a }} = 2 9w w ^ 9a9p ~ 9a9(s > ' 

V w-» = I} g™ cos 9 W (g^g - g^g J ) , V w - a0 = ^9w cos 6 W {9 a 90 - 9^9 fi ) ■ 



(26) 



Corresponding vertices exist for W + / W~ bosons. The decomposition of four particle vertices involving 
the Higgs boson introduces a new scalar pseudoparticle, which we denote by h^. In order not to generate 
fake four particle vertices we define it not to be self-conjugate. The corresponding vertices read 

■,*h,h,h -»,h,/i4 p -,,h.h 
V h ~ * V h ' ^4 ' v hi ' 
■\)h, Z n, Zv -,,h,hi p -,,Zp,,Zv f07\ 
v h ~ * v h ' r hi ■ Yht ' y Z< ) 

where the interactions of the /14 pseudoparticle are defined by 

Vh, h 4 
h =*. 

^4 _ * w 2 ' 

Z^,Z„ = ■ gg, .... ( 28 ) 
hi 1 2 cos 2 e w 



y hi - i 2cos ^ 9 

+ n 2 

^ - 11V 

v hi — 1 2 9 ' 

and where we have introduced the scalar "propagator" of the /14 pseudoparticle 
P hi = i ■ (29) 

Since all remaining vertices in the Standard Model are three point vertices, the vertex decomposition is 
hereby complete. The additional Standard Model propagators and vertices arising from this decomposition 
are summarised in Tabs. [Hand respectively. 
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g,e,KL 



g,e',MN 



94, HG 
h h 



V2 



VVT (£,£') 



hi 



V 




W~,E 



9 X 2 



VVS (e, e') where 



= i g w VVT (e, e') 



Xw = 1 
Xz = cos 6w 



W~,E 



A/Z, e' 



i g w k-a/z VVT (e, e') where 



= siu&vc 



WS (e, e') = , VVT^" (e, e') = | (,g^ A 5 ^ - ^".g^) e A e' K 



Tab. 2 Standard Model vertices arising from the vertex decomposition and replacing the 
four particle vertices. In this context e and e' denote arbitrary incoming vector 
currents. Note that due to the antisymmetry of VVT' 1 ", we can make the replace- 
ment D^WT^ = 2VVT M ", which leads to a slight decrease in computation 
time. 



2.5 Prefactors of diagrams with external fermions 

When calculating currents with an arbitrary number of possibly indistinguishable external fermions, we have 
to take into account, that each Feynman diagram contains a prefactor 



5 = (-l) s ^ CT1 ' 



(30) 



according to the number of fermion permutations S/ in the external particle assignment a — (ax, ...,a n ). 
To be used in the context of a recursive computation, this prefactor must be defined on a local basis in order 
to avoid the proliferation of information on different a. It is then sufficient to note that Eq. (|30p holds at 
the level of interaction vertices. More precisely we can define the local prefactor S (tti, 7^) of Eq. © as 

5(7ri,7r 2 ) - (-l) 8 /^'"*) . (31) 

Here S/ (7^,772) counts the number of fermion permutations that is needed to restore a predefined, for 
example ascending index ordering when combining the sets tti and 7T2 into the set 7r = 7Ti ® 112- Upon 
iterating this procedure, we obtain the correct relative prefactors S for each diagram. 



3 Matrix element generation in Comix 

The general formulae to recursively compute tree-level amplitudes have been stated in Sec.[H Here we briefly 
explain, which conventions are used to define the external particle currents and internal Lorentz structures. 
We also elaborate on how to organise the computation and how to reduce the effective computation time 
per phase space point by a multi-threaded structure of the implementation. 
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Fig. 1 Structure of the multi-threaded implementation for matrix element com- 
putation in COMIX. The number of threads N is variable and depends on 
the number of available processors. The main program communicates start 
and wait signals to the calculator threads, while those communicate done 
and wait signals to the main program. Details are explained in the text. 



3.1 Choice of the spinor basis 

We employ the spinor basis introduced in Ref. [23] . Accordingly, the 7-matrices are taken in the Weyl 
representation. The main advantage of this representation is that spinors for massless particles are described 
through two nonzero components only. This fact greatly alleviates their construction as well as the evaluation 
of vertices. Polarisation vectors for external vector bosons are constructed according to Ref. [21]. As pointed 
out in Sec. [2j within the Standard Model tensor particles never occur as external states, such that there is 
no need to explicitly construct polarisation tensors. 



3.2 Implementation details 

The algorithms presented in this paper are intended to be used for large multiplicity matrix element calcula- 
tions. In this context, it is often useful to sample over helicities of external particles in a Monte Carlo fashion. 
However, this introduces additional degrees of freedom and leads to a slower convergence of the integral. 
Furthermore when taking Eq. © serious, we note that for helicity-summed ME's, it is possible to reuse 
currents to compute amplitudes with different configurations. Namely if the helicities of external particles 
assigned to a particular current do not change, it does not need to be recomputed. This leads to a significant 
decrease in evaluation time for the helicity summed ME's compared to the naive method of computing the 
full amplitude afresh for different configurations. A corresponding comparison can be found in Sec. El The 
default choice in COMIX is helicity summation. To allow computations for very large multiplicities, however, 
helicity sampling can be enabled as an option. 

The effective computation time per phase space point can be further reduced by a multi-threaded implemen- 
tation of Eq. |6]). Figure Q] shows the basic structure of this algorithm. The main advantage of Eq. |6|) is, that 
in order to compute a current that depends on n external particles, it is sufficient to know all subcurrents 
that depend on to < n external particles. This leads to a straightforward multi-threading algorithm. 

• Create N threads at program startup with the following properties 

1. The thread waits for the main program to signal the start of a computation. 
It then signals the main program to wait. 

2. It takes a number n and computes a block of currents depending on n external particles using 
subcurrents depending on m < n external particles. If n = 1, it computes external polarisation 
vectors and spinors. 

3. It signals the main program that the calculation is done and returns to step[TJ 
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• For each phase space point, employ the following algorithm in the main program 

1. Start with n = 1. 

2. Split the number of currents that depend on n external particles into N blocks. 
Communicate n and one block to each calculator thread. 

3. Signal the threads to start their computation. 
Wait for all threads to signal completion. 

4. Let n — > n + 1 and return to step [2] if further currents need to be computed. 

The efficiency of this algorithm solely depends on an efficient thread library. The overhead with a modern 
POSIX threading is about 10% of the total computational cost. This, however, is not of any concern 
considering that the employment of multiple CPU's reduces the computation time roughly proportional to 
the increase in processor usage. 



4 Integration techniques in Comix 

In this section we present two new methods for integrating over the phase space. Both of them are designed 
to cope especially with large numbers of outgoing particles. The first method is a fully general approach and 
makes use of the standard multi-channel technique |26] in a recursive fashion, i.e. the phase space sampling 
fits the method of generating the corresponding matrix element. The second method is designed for QCD 
and QCD-associated processes and employs the phase space generator HAAG [19] in conjunction with a new 
prescription for coupling colour and momentum sampling and the multi-channel technique. 

4.1 Recursive algorithm for phase space integration 

One of the most efficient general approaches to sample the phase space of multi-particle processes is, to 
employ a multi-channel method according to Ref. |26j with each of the single channels corresponding to 
the pole structure of a certain Feynman diagram. However, for large numbers of diagrams this is clearly 
not the method of choice. In the following we will therefore focus on the recursive relations for phase space 
generation proposed in Ref. |18j . We construct a separate multi-channel for each possible subamplitude on the 
flight according to the propagator structure and use VEGAS |27j to optimise the integration over propagator 
masses and polar angles in decays. The obvious drawback of this procedure is evident: It relies heavily on 
the assumption that the matrix element factorises according to its propagator structure. However, it is a 
generalisable way to tame the rather factorial growth in the number of phase space channels encountered in 
conventional approaches [8,9, 10 . If we take the prescription serious, we can factorise the full phase space 
weight such that it can be computed in a recursive fashion corresponding to how the matrix element is 
evaluated. It turns out that this procedure gives an excellent performance, cf. Sec. 03 



4.1.1 Brief review of phase space factorisation 



In the following we consider a 2 — > n scattering process and denote incoming particles by a and b and 
outgoing particles by 1 . . . n. The corresponding n-particle differential phase space element reads 



d<J>„ (a, 6; 1, . . . ,n) 



n 



\ (2?r) £ 



5 (p* - mf) 9 (pa,) 



(32) 



where are the on-shell masses of outgoing particles. Following Ref. | 
factorised according to 

d$ n (a, b; 1, . . . ,n) = d$„_ m+ i (a, b; %, m + 1, . . . , n) d$ m (%; 1, . . . , m) , 

Z7T 



the full phase space may be 



(33) 



where it = {1, . . . , m} corresponds to a set of particle indices, similar to Sec. [H Generally in this section, 
Greek letters denote a subset of all possible indices. Overlined letters denote the missing subset, i.e. a = 
{a, b, 1, . . . , n} \ a for all a C {a, b, 1, . . . , n}. Similarly, W1W2 = ( {a, b, 1, . . . , n} \ tt%) \ tt2, etc. Equation |33|) 
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Fig. 2 Basic vertices for phase space generation. Grey blobs correspond to eventually off mass-shell 
particles. Dark blobs denote known momenta, light blobs unknown momenta. Arrows indicate 
the momentum flow, i.e. the order in which unknown momenta are determined from known ones. 
The Z)-vertex corresponds to overall momentum conservation. 



allows to decompose the complete phase space into building blocks corresponding to the t- and s-channel 
decay processes T* b ab7T = d<I>2 (ce, b; it, abir) and S£'*^ p = d<£>2 (t; p,n\ p) and the s-channel production 
process D a ^ , cf. Fig. [21 We refer to these objects as phase space vertices, while the integral P K — ds- R /2 / K, 
introduced in Eq. ([SU]) . will be called a phase space propagator. We use the same notation as for the 
propagators in Seed] to highlight the close correspondence between matrix element computation and phase 
space generation. In the algorithm presented here, only timelike propagators are employed. 

The phase space vertices are used differently in the case of weight calculation and phase space generation. 
Consider the i-channel decay. If a phase space point is to be generated, the new final state momenta p v and 
p^j^r are determined from the known initial state momenta p a and p b - If a weight needs to be computed, the 
new weight Wa is determined from the vertex weight and the input weights and w-^^. The corresponding 
situations are depicted in Figs. [5] and [31 respectively. The basic building blocks of phase space integration 
are summarised as follows 

{1 if 7T or 7r external 
dsjr , 
else 
2tt 

^-^^W, (34, 

£>„,, = (2») 4 d'p- t S»> (p„ + p„ - pjs) . 
Here we have introduced the triangular function 



A (s a , s b , s c ) = \j (s a - s b - s c ) 2 - 4s b s c (35) 

Note that even since a might correspond to an off-shell internal particle, b always indicates a fixed external 
incoming particle. This is essential in all further considerations and allows reusing weight factors in the 
Monte Carlo integration, just as currents are reused in the matrix element computation. The functions 
corresponding to and T^' abK are in fact identical, since they represent a solid angle integration. In 

practice however we choose the different sampling strategics proposed in Rcf. 18J . The s-channel production 
vertex D a b is only needed for bookkeeping, since it corresponds to overall momentum conservation and the 
associated overall weight factor (27r) 4 . 



4.1.2 A simple example 

We illustrate in this section how a recursive phase space generator for the process qq —> e + e~g can be 
constructed, based on the diagrammatic structure of the integrand. Figure H] depicts the translation of the 
corresponding Feynman diagrams into related building blocks of the phase space. The standard procedure 
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Fig. 3 Basic decay vertices for weight calculation. Dark blobs denote potentially nontrivial known 
weights, light blobs weights to be determined. Arrows indicate the weight flow, i.e. the order in 
which unknown weights are determined from known ones. The D-vertex corresponds to overall 
momentum conservation. 



to define an integrator consists of constructing one integration channel per line of Fig. 2] and joining these 
channels in a multi-channel. Because it is based on full diagrams, this strategy cannot be implemented in a 
recursive fashion and we have to modify it. Consider first, which tasks have to be performed for each phase 
space point. 

To generate momenta, one starts with the s-channel propagator P23- Then, depending on what yields the 
better performance, the t-channel decay T L or T b and finally the s-channel decay S^ 3 are constructed, 
leading to the final state momenta p\ . . .p%. Note again that D-type vertices are just for bookkeeping at 
this point. They simply imply overall momentum conservation. When computing the phase space weight, 
the order of treating vertices can be reversed because all corresponding momenta are known. Therefore the 
weights for the decay S^ 3 and the propagator P23 are computed first, followed by the weig hts for T^f 3 and 



T ab ' . It is obvious that the weights P23 and S 2 3° are unique and therefore have to be computed only once, 
although arising in both lines of Fig. HI We refer to this feature as the "weight flow" , which is directed from 
the final state particles and the right beam, particle b, towards the left beam, particle a. This generalises to 
arbitrary processes, provided that the right beam particle is kept fix, i.e. i-channel indices always combine 
only a and external indices, as indicated in Figs. [5] and [31 It allows to compute the full phase space weight 
recursively, in a manner similar to Eq. ([H]), which implies in particular, that at most the same growth is 
induced in the matrix element and the phase space weight computation. 

Let us illustrate this new procedure using the above example. The difference with respect to the standard 
approach is how multi-channels are defined. Following the weight flow, in the first step of the recursion we 
construct a multi-channel for the phase space element d<I>2 ({23}; 2, 3). Of course, since particles 2 and 3 are 
external, this multi-channel consists of one single channel only, which is the s-channel decay S^ 3 - It has 
therefore no additional parameters. In the second step we construct a multi-channel for the full phase space 

X 23 23 1 

d$3 (a, b; 1,2, 3), which receives contributions from the two t-channels T a ' b and T a b ' . Each of them can be 
assigned a multi-channel weight w, which eventually yields the overall weight 



,2,3 



(2tt) 4 F- 



1,23 T7. 

W a[b F 


rji 1,23 i ) 2,3 
I a,b ^23^23 


1 23,1 Tji 
+ W a,b F 


rji 23,1 / > 2,3 
I a,b -^23^23 




1.23 

W a,b - 


23,1 
^ W a,b 





(36) 



Here P denotes the propagator weight, and F is a generalised mean function, see next section. The overall 
factor (27r) 4 arises from the -D-type vertices. 

We are left with the task to determine the necessary building blocks of the phase space using information 
from the matrix element. This turns out to be extremely simple since we now in fact have a phase space 
weight recursion of the form of Eq. ([6]) (cf. Eqs. (|38f and ([39]) ). where multi-channels are associated with 
intermediate s/t-channcl propagators. Therefore each intermediate current in the matrix element implies 
a separate multi-channel in the corresponding integrator and each vertex implies a decay vertex associated 
with a single channel. Because of this correspondence the default sampling strategy for s-channel propagator 
masses can be chosen according to the type of intermediate particle, i.e. a Breit-Wigner-like distribution for 
massive, unstable particles and a l/s a -type distribution for massless particles. These distributions as well as 
the polar angle distributions in decays (cf. Eq. fSj]) ) are further refined during the integration using VEGAS. 
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4.1.3 Formulation of the recursive algorithm 



In this subsection we derive the general algorithm for the recursive phase space integrator. We employ the 
notation of Sec. 14.1.11 Recursive relations for phase space integration in terms of the quantities introduced 
ibidem can be defined through 



d$ s (tt) = P W1 d$ s (tti) P W2 d$ s (7r a ) 

d4 b) (a) = TlX** P ni d$ s (tti) P^ (otti) 



(iri,ir 2 )eOT'(ir) 



+ A*,b d$ S (55) 



(37) 



The above equations correspond to selecting one possible splitting of the multi-index n or ab per phase space 
point. We can improve the integration procedure by forming an average over all possible splittings in the 
spirit of a multi-channel. Let F be a generalised mean function. We can then use the F-mean to define 



d$ s (tt) = F- 1 



E < 1,7r2 

v (ir I ,ir a )eOT'(jr) 



x £ w* 1 '^ F [ S?>** P 7Tl d$ s (tti) d$ s (7r a ) 

(7ri,7r 2 )ee)P(7r) 



(38) 



d^ (a) = F- 1 



+ E 



F 



(7ri,7r 2 )ee>'P(c(6) 



D Qi6 d$ s (a6) 



E 



v 



x F 



T^ 2 P W1 d$ s (Tn) P X2 d$^ fc) (cwri) 



(39) 
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In this context we define the one- and no-particle phase space 



d$ (i) = l, 
d$ (0) = . 



(40) 



The function u> corresponds to a vertex-specific weight which may be adapted to optimise the integration 
procedure, see Ref. |26j . The second sums run over all possible S- and T-type vertices which have a corre- 
spondence in the matrix element. The full differential phase space element is given by 



Note that Eqs. (|3"5|) and (|3"§|) in the form stated above are not suited to generate the sequence of final state 
momenta. To do so one rather has to employ the following algorithm, which corresponds to a reversion of 
the recursion and respects the weight factors w introduced above. 

• From the set of possible vertices connecting currents in the matrix element, choose a sequence con- 
necting all external particles in the following way: 

1. Start with the set of indices n = {b, 1, . . . , n}, 
corresponding to the unique external current of index a. 

2. From the set of possible phase space vertices connecting to it select one according to an on the flight 
constructed multi-channel employing the weights w0 If ir is a single index, stop the recursion. 

3. According to the selected vertex, split ir into the subsets m and tt2- Repeat step [2] for these 
subsets. 

• Fore each vertex, make use of the fact that it is equivalent to W and adjust the indices in an appropriate 
way for momentum generation. That is if any n contains b and other indices, replace tt by tt. 

• Order T-type vertices ascending and 5-type vertices descending in the number of external indices 
connected to initial states. 

• Generate the corresponding momenta starting with T-type vertices. 

Even though T-type vertices depend on 6, since b is fixed throughout the computation of one phase space 
point we obtain no expressions depending on more than two particle indices. This induces the same growth 
of computational complexity in both the hard matrix elements and the phase space and makes the above 
algorithm well suited for integration of processes with large final state multiplicity. In the following we refer 
to it as the Recursive Phase space Generator (RPG). 

4.1.4 Implementation details 

Since the phase space weight computation, Eq. (|38|) obeys a recursion similar to those of the matrix element 
calculation, Eq. ([5]), it is straightforward to implement this weight computation into a numerical program 
along the lines of Sec. 13.21 The same techniques described for the multi-threading of matrix element calcu- 
lations can be implemented for the phase space weight. In the multi-threaded version of COMIX, this weight 
is computed in parallel to the matrix element, which further reduces the net computation time if enough 
resources are available. 

4.2 Colour sampling 

For QCD and QCD associated processes with a large number of external legs, it becomes unfeasible to 
compute colour-summed scattering amplitudes. Instead the better strategy is to sample over external colour 
assignments in a given representation of SU(3). According to Eqs. (O - (fTO]) . this selects a set of colour- 
ordered amplitudes which contribute to the corresponding point in colour space. This set is typically strongly 
reduced compared to the full set of partial amplitudes. The issue has been studied in Ref. [55] for the 
fundamental representation decomposition, the adjoint representation decomposition and the colour flow 

2 Note that in this context weights have to be normalised to unity on the flight. 



d$ n (a, 6; 1, . . . ,n) = d$ T (a) . 



(41) 
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decomposition, which has been presented therein. The conclusion is that the colour flow decomposition is 
the method best suited for sampling over colour assignments if the number of external partons is large, i.e. it 
provides the slowest growth in the average number of partial amplitudes per non- vanishing colour assignment. 
Also it has been exemplified for recursive calculations in Ref. [15], that the colour flow decomposition is 
advantageous, since no computational intensive matrix multiplications have to be performed. We therefore 
employ this prescription throughout COMIX. 

In the following we focus on an n-gluon scattering process. However, the presented ideas and algorithms are 
straightforward to generalise for arbitrary sets of colour octet objects, such as e.g. quark-antiquark pairs. In 
the colour flow decomposition each external gluon is labeled by a colour index i and an anti-colour index j. 
The colour assignment for an n-gluon scattering is thus given by selecting each index ii, . . .i n and Ji, . . . J„ 
out of three values (R, G, B) and (R, G, B) . 

4.2.1 Determination of colour flows from colour assignments 

A specific colour flow, and thus an ordering in the sense of a colour-ordered amplitude, is specified by a 
permutation 

a=(l,a- 2 ,a 3 ,...,(X n ) € S n -i (42) 
of external gluon indices. This colour flow contributes to a colour assignment, if 

§ilJ*1 §i<r 2 3<r 3 . . . £V„Jl = 1 _ (43) 

It is thus easy to construct an algorithm which determines all valid colour flows from a given colour assign- 
ment. 

1. Set the first gluon index to o\ — 1. Let k — 2. 

2. Select one of the remaining gluon indices to be o^, such that i - fc _ 1 = j ak . If this is possible, let 
k — ► k+ 1. Otherwise let k — ► k — 1, then repeat this step selecting a different <7fc. 

3. If k = n + 1 and i an = j ai , a valid flow has been found. 
Otherwise continue with step [U 

By systematically selecting through all possible at in step [2] all valid colour flows are determined. 



4.2.2 Selection of colour assignments 

The simplest way of choosing a colour assignment is accomplished by randomly selecting the 2n colours for 
the i- and J- indices. Each colour is chosen with an equal probability, leading to a weight of 3 2n . However, 
only a small fraction of those assignments will have at least one valid colour flow. A trivial (but not sufficient) 
condition for non-vanishing amplitudes is, that the number of i-indices carrying the colour R (G,B) must 
be equal to the number of j-indices carrying the corresponding anticolour. 

We thus propose a more efficient way to determine colour configurations. 

1. The n «-indices are selected randomly in (R, G, B). 

2. A permutation a — (a±, . . . , a n ) of n particles is selected randomly with a uniform weight. The 
anticolours of the J- indices are then given by jk = V fc , for fc = 1, . . . , n. 

3. Each colour assignment is weighted by 

w = 3" , U \ , , (44) 
n R \n G \n B \ 

where n^, uq and tib are the multiplicities of i-indices 
carrying the colours R, G and B, respectively. 

Clearly, assignments generated by this algorithm will always fulfil the trivial condition mentioned above. 
Moreover, the weight is roughly proportional to the number of possible colour flows and thus already corre- 
sponds to some extent to the expected cross section for this colour configuration. 
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4.2.3 A simple example 



To illustrate the colour sampling in the colour flow decomposition we consider a five gluon scattering process. 
The random selection of a colour configuration using the improved algorithm may return the following i- 
indices: 

h = R, i2 = R, h = G, u = G, i 5 = B. (45) 
The j-indices are fixed by a randomly chosen permutation, say a — (4, 1, 2, 5, 3): 

Ji = G, j 2 = R, 33 = R, H = B, 3i=G. (46) 
For this assignment the only colour flow that satisfies Eq. ([4*31 is given by the permutation a = (1, 4, 5, 3, 2). 

4.3 Combined colour-momentum integration techniques 

Generally the peaking behaviour of the colour-sampled differential cross section is rather complex within 
the phase space and strongly different for different colour assignments. The idea must thus be to construct 
integrators specific for a given colour assignment, based on the knowledge of contributing partial amplitudes. 
One can for example think of a variant of the algorithm described in Sec. 14. 11 where the basic building blocks 
of the phase space are either available or not, depending whether there is a corresponding non- vanishing 
coloured current present in the matrix element. However, in practice this choice does not lead to any 
significant improvement of the integration behaviour of the RPG and we thus refrain from promoting this 
method^ Instead we present a second type of integrator, dedicated to be used with QCD and QCD associated 
processes, which is based on the HAAG algorithm 19J. As before we concentrate on purely gluonic processes. 



4.3.1 Integration of partial amplitudes and colour configurations 

As a basic building block we use the HAAG-integrator, which generates momenta distributed proportional 
to a QCD antenna function [19], 

A n (p!,p 2 , ...,p n ) = — : -. (47) 

\PlP2 ) (P2P3 )-{Pn-lPn) (PnPl ) 

Details on our implementation of the algorithm and improvements to the original version are given in Ref. [7] . 
Single HAAG-channels provide efficient integrators for squared partial amplitudes associated with a given 
colour flow, both labeled by the same permutation a, Eq. (|42l) . For the HAAG-channel the permutation 
corresponds to the order of momenta in the antenna function. As for the RPG we again obtain a close 
correspondence between the matrix element and the phase space generation, now at the level of partial 
amplitudes. 

The cross section for a single colour assignment is given by the squared sum of partial amplitudes associated 
with all valid colour flows. Ignoring the interferences between the amplitudes in the context of the phase space 
setup, a dedicated integrator can be constructed by combing the corresponding HAAG-channels for each flow 
in a multi-channel integrator. With growing number of external particles, however, one faces the following 
problem: Although the average number of contributing colour flows per colour assignment is relatively low 
in the colour flow decomposition, the maximal number grows factorially. Thus it quickly becomes impossible 
to store all data associated with the multi-channels defined above, i.e. the contributing HAAG-channels and 
the internal weights. The situation gets even worse if the sampling over all colour assignments is taken 
into account, because the number of possible assignments grows exponentially with the number of external 
particles. The solution to this is not to store any multi-channel parameters, but to generate the complete 
multi-channel on the flight. 

A fast algorithm, as presented in Sec. 14. 2. II to provide all colour flows from a colour assignment is essential 
for this approach: for a single phase space point one has to loop three times over the list of all colour flows 
(which cannot be stored as well due to the factorially growing maximal number of flows). 

1. To determine the normalisation of weights within the multi-channel integrator. 



3 Note that this is not a statement about the integration behaviour of the RPG itself, but only about a possible coupling of 
colour and momentum sampling using the RPG. 
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2. To select a channel for generating a phase space point with a probability given by the relative weight 
atk, and 



3. To compute the multi-channel weight corresponding to this phase space point. 

4.3.2 Optimisation techniques 

The proposed integrator contains a number of parameters which can be adjusted or adapted to reduce the 
variance during integration. A multi-channel integrator dedicated to a specific colour assignment has the 
following degrees of freedom for optimisation: 

• VEGAS grids to refine individual HAAG-channcls, 

• Relative weights ctk in the multi-channel generator, 

The sheer multiplicity of different channels and on-the-flight construction of the integrator forbids an indi- 
vidual adaptation of all parameters. However, their number can be greatly reduced by making use of the 
symmetry among different HAAG-channels w.r.t. to permutations of the final state. All channels with the 
same relative positions of the initial state momenta within the antenna can be determined from each other 
by a permutation of final state momenta. This prevents the number of structurally different phase space 
channels from growing factorially with the number of particles and induces a linear growth only. Taking 
into account that the same symmetry holds for the partial amplitudes justifies to reuse the optimisation 
parameters among all channels of one kind. For later reference we label different types of HAAG-channels 
(and respective partial amplitudes) by the number of final state momenta between the first and the second 
incoming momentum within a certain antenna. 

We achieve the best integration efficiency by performing the optimisation of the free parameters prior to 
the actual integration: The VEGAS grids of the HAAG-channels are adapted individually by integrating 
corresponding single squared partial amplitudes over the allowedphase space. Using the above mentioned 
symmetry this has to be done only for one channel of each kincQ This technique not only speeds up the 
optimisation, it also provides a much cleaner environment for the adaptation of the VEGAS grids. At this 
stage a summation over helicities is performed. Cross sections at, given by the integration result from the 
channel of type t, are stored. 

In the actual integration run no further optimisation is performed. The channels are used as they emerged 
from the above procedure, including the VEGAS-grid and a parameter a^, proportional to the cross section, 
a tl of the corresponding squared partial amplitude. 

Best performance is achieved, if the colour assignment is selected with a probability proportional to the 
sum of cross sections of contributing squared partial amplitudes (as determined during the optimisation 
step), instead of the weight given by Eq. (JUJ). To do so, the total normalisation for the new weight must be 
determined summing over all colour assignments. For n-gluon processes this number is given by the following 
simple formula: 



where the (Tjnia.(i,n—i—2) is the cross section of a squared partial amplitude of the type "min(i,n — i — 2)". 
The reweighting can be done by a simple hit-or-miss method. 

For the integration run it is a matter of choice whether to sum or sample over helicities. All practical tests 
for up to the 11-gluon process favoured summation. Beyond that, however, it seems to become too costly to 
compute summed matrix elements, thus a sampling should be considered. 

In the context of this work, we refer to the above algorithm as the Colour Sampling Integrator (CSI). 



4 During this step the full result can not be determined since potential interferences between partial amplitudes are ignored. 
However, it is sufficient for computing the leading X/Nc limit for n gluon processes, using the fact that in the colour flow 
decomposition (as well as in the fundamental representation decomposition) interferences are always sublcading. 



n-2 




(48) 



1G 



5 Results 



In this section we present selected results generated with COMIX. We focus on the special feature of this 
new generator, to be suitable in particular for computation of large multiplicity matrix elements. A detailed 
comparison of integration times, compared to a dedicated code using CSW vertex rules and the generator 
AMEGIC++ can be found in Ref. [7j. 

5.1 Helicity summation vs. helicity sampling 

Firstly we illustrate the effect of suitable matrix element generation in the helicity summed mode of COMIX, 
see Sec. 13.21 Computation times for helicity summed and helicity sampled matrix elements in pure gluonic 
processes are compared in Tab. [3l The naive ratio between the two is the number of possible helicity 
assignments of the respective amplitude, 2™ — 2(n + 1), with n the number of external gluons. This naive 
ratio corresponds to computing the amplitude afresh for each of the different helicity assignments. Employing 
the ideas presented in Sec. 13. 2\ however we find that this value overestimates the real computational cost by 
up to a factor of 7. Obviously this statement is process dependent. The general feature, however is that 
there is a gain when computing helicity summed matrix elements. For the computation of cross sections 
this type of calculation might be preferred over the helicity sampled mode, especially when using the phase 
space integration methods of the previous chapter, which are not designed for helicity sampling. 

5.2 Performance of the CSI and 2 — > n gluon benchmarks 

In this subsection we present a comparison of gluon production cross sections to illustrate both the perfor- 
mance of the CSI and the efficiency of the matrix element generation. We start with a fixed centre-of-mass 
energy. The parameters are those of Refs. [19,22 , i.e. as = 0.12 and 

p Tl >60GeV, M<2, ARij>0.7, (49) 

for all final state gluons i and pairs of gluons Integration results are summarised in Tab. [H We find 
perfect agreement with the results in the literature and give new predictions for the processes gg — > llg and 
gg — ► 12g. Results have been generated with the CSI, except for the 2 — > 11 and 2 — * 12 process, where 
RAMBO 30J has been employed. In order to examine the performance of the new phase space generator in 
a more realistic scenario, we investigate the same partonic processes at the LHC and employ the Tevatron 
Run II \lt algorithm [3ljl to define a cut on the multi-particle phase space. The respective results are 
summarised in Tab. We find that the CSI performs very well in both cases, even for large multiplicities, 
such that the respective cross sections can be computed with good precision. 

Figures [5] and [S] show the convergence behaviour of the CSI for various gluon multiplicities. Since the 
computation of 2 — * 8 and 2^9 gluon processes is quite cumbersome, it is worthwhile to switch to the 
helicity sampled mode in that case. Correspondingly we compare the performance of the CSI in helicity 
summed and helicity sampled mode in Fig. [6] 

5.3 Performance of the RPG and comparison with other generators 

We finally compare the performance of COMIX with those of other programs. All results presented in 
this section were obtained with the RPG described in Sec. 14.11 As references we use AMEGIC++ [5] and 
ALPGEN pj]. The original setup for the comparison was established during the MC4LHC workshop 32J. 
For a comprehensive listing of results from all participating projects, see ibidem. Input parameters are given 
in Tab. [51 Cross sections are summarised in Tabs. 171-191 and ITD1 

As pointed out in Sec. 14.11 a drawback of the RPG is that it might not be able to adapt to certain peaks of 
the matrix element which correspond to specific diagrams. No significant disadvantage compared to other 
generators can however be observed. A measure for the efficiency of a phase space generator is given by 
the ratio of the average over the maximal weight (w)/u> max , i.e. the efficiency for generating events of unit 
weight using a hit-or-miss method. As discussed in Ref. [33j . the maximum weight and thus this ratio is 
a numerically rather unstable quantity, often determined by very rare events in the high tail of the weight 

5 Note that we replace Ai?? — » 2 [cosh Ajjij — cos A<f>ij ] in order to match the Durham measure for final state clusterings. 
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distribution. In Tab. 1111 we therefore list the more stable quantity (u>)/u>^ ax , where the reduced maximum 
weight w^ ax is defined such that 1 — (min(u>, w^ nax ))/('w) = £«1. It turns out that we achieve a reasonably 
good efficiency using the RPG, even for very large multiplicities. It can therefore be concluded that this 
phase space generator is an excellent approach to tame the factorial growth of phase space channels while 
still maintaining an a priori adaptation to the assumed peak structure of the integrand. 

6 Conclusions 

We have presented the new matrix element generator COMIX, based on the colour dressed Berends-Giele 
recursive relations and two new methods for phase space generation. We have analysed the performance 
of the new generator and compared the respective results to other ME generators. We find that the new 
algorithms perform very well and we obtain promising results for large multiplicity processes. COMIX can 
therefore be considered an excellent supplementary generator for large multiplicities, which is especially 
helpful in the context of a matrix element - parton shower merging. The treatment of colour in COMIX 
makes the algorithm well suited for such an interface, since the colour structure of the matrix element does 
not need to be guessed from the kinematics, it is rather fixed on a point by point basis. A corresponding 
publication is forthcoming [34J. 
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Tab. 3 Computation time for multi-gluon scattering matrix elements sampled over colour con- 
figurations. Displayed times are averages for a single evaluation of the colour dressed BG 
recursion relation, when summing and sampling over helicity configurations, respectively. 
Additionally in the last column, labeled 'Gain' we give the inverse ratio of evaluation 
times multiplied by the naive ratio 2™ — 2(n + 1), where n is the number of external 
gluons. Numbers were generated on a 2.80 GHz Pentium® 4 CPU. 
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Tab. 4 Cross sections for multi-gluon scattering at the centre-of-mass energy ^/s, using the 
phase space cuts specified in Eq. (149 1 , compared to literature results. In parentheses the 
statistical error is stated in units of the last digit of the cross section. 
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Comix 


2703(14) 


407.0(36) 


66.5(13) 


15.2(26) 



Tab. 5 Multi-gluon cross sections at the LHC with yd > 20 GeV and d defined as in Ref. [31] . 

except that Aiiy —* 2 [ cosh Arjij — cos A(/>y ] . In parentheses the statistical error is 
stated in units of the last digit of the cross section. 
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Fig. 5 Overall integration performance for multi-gluon scattering. Upper panels display the Monte Carlo 
estimate of the cross section with the corresponding la statistical error band as a function of the 
total integration time. Lower panels show the relative statistical error. HAAG denotes the phase 
space integrator described in Ref. [7], applied on colour- and helicity-summed ME, generated using 
the CSW vertex rules. CSI denotes the integrator discussed in section 14.3.11 applied on colour- 
sampled and helicity-summed ME's, generated using the CDBG recursion. Results for RAMBO 
were generated using colour- and helicity-sampled ME's form the CDBG recursion. Calculations 

TM 

have been performed on a 2.66 GHz Xeon CPU 
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Fig. 6 Overall integration performance for multi-gluon scattering, continued from Fig. [5] Additionally, for 
the CSI a sampling over helicity is considered, denoted by CSI(HS). 



Parameter 


Value 


EW parameters in the scheme 


G F 


1.16639 x 10" 5 


OiQED 


1/132.51 


sin 2 6w 


0.2222 


M w 


80.419 GeV 


M z 


91.188 GeV 


m H 


120 GeV 


CKM matrix 


V u d, v cs 


0.975 


QCD parameters 


PDF set 


CTEQ6L1 




0.130 




M z 


jet, initial parton 


g, u, d, s, c 



Parameter 


Value 


Non-zero fermion masses (no evolution) 


m b 


4.7 GeV 


m t 


174.3 GeV 


m T 


1.777 GeV 


Widths (fixed width scheme) 




2.048 GeV 


T z 


2.446 GeV 


r H 


3.7 x 10~ 3 GeV 




1.508 GeV 


T T 


2.36 x 10~ 12 GeV 


Cuts 


P±,i 


> 20 GeV 


\m\ 


< 2.5 




> 0.4 


no cuts on particles of to > 3 GeV and Vi 



Tab. 6 Parameters for the MC4LHC comparison setup. 
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a [fib] 


Number of jets 


jets 


2 


3 


4 


5 


6 


7 


8 


Comix 


331.0(4) 


22.72(6) 


4.95(2) 


1.232(4) 


0.352(1) 


0.1133(5) 


0.0369(3) 


ALPGEN 


331.7(3) 


22.49(7) 


4.81(1) 


1.176(9) 


0.330(1) 






AMEGIC 


331.0(4) 


22.78(6) 


4.98(1) 


1.238(4) 









a [/ib] 


Number of jets 


bb + jets 





1 


2 


3 


4 


5 


6 


Comix 


471.2(5) 


8.83(2) 


1.813(8) 


0.459(2) 


0.150(1) 


0.0531(5) 


0.0205(4) 


ALPGEN 


470.6(6) 


8.83(1) 


1.822(9) 


0.459(2) 


0.150(2) 


0.053(1) 


0.0215(8) 


AMEGIC 


470.3(4) 


8.84(2) 


1.817(6) 











a [pb] 


Number of jets 


tt + jets 





1 


2 


3 


4 


5 


6 


Comix 


754.8(8) 


745(1) 


518(1) 


309.8(8) 


170.4(7) 


89.2(4) 


44.4(4) 


ALPGEN 


755.4(8) 


748(2) 


518(2) 


310.9(8) 


170.9(5) 


87.6(3) 


45.1(8) 


AMEGIC 


754.4(3) 


747(1) 


520(1) 











Tab. 7 Cross sections a in the MC4LHC comparison [32] setup. In parentheses the statistical error is stated 
in units of the last digit of the cross section. Note that for AMEGIC++ and COMIX all subprocesses are 
considered, while ALPGEN is restricted to up to four quarks. 



a [pb] 


Number of jets 


e + v e + jets 





1 


2 


3 


4 


5 


6 


Comix 


5434(5) 


1274(2) 


465(1) 


183.0(6) 


77.5(3) 


33.8(1) 


14.7(1) 


ALPGEN 


5423(9) 


1291(13) 


465(2) 


182.8(8) 


75.7(8) 


32.5(2) 


13.9(2) 


AMEGIC 


5432(5) 


1279(2) 


466(2) 


185.2(5) 


77.3(4) 







a [pb] 


Number of jets 


e~v e + jets 





1 


2 


3 


4 


5 


6 


Comix 


3911(4) 


1011(2) 


362(1) 


137.1(3) 


54.9(2) 


22.4(1) 


9.26(4) 


ALPGEN 


3904(6) 


1013(2) 


364(2) 


136(1) 


53.6(6) 


21.6(2) 


8.7(1) 


AMEGIC 


3903(4) 


1012(2) 


363(1) 


137.6(3) 


54.8(6) 







a [pb] 


Number of jets 


e~e + + jets 





1 


2 


3 


4 


5 


6 


Comix 


723.5(4) 


187.9(3) 


69.7(2) 


27.14(7) 


11.09(4) 


4.68(2) 


2.02(2) 


ALPGEN 


723.4(9) 


188.3(3) 


69.9(3) 


27.2(1) 


10.95(5) 


4.6(1) 


1.85(1) 


AMEGIC 


723.0(8) 


188.2(3) 


69.6(2) 


27.21(6) 


11.1(1) 







a [pb] 


Number of jets 


VeP e + jets 





1 


2 


3 


4 


5 


6 


Comix 


3266(3) 


715.9(8) 


266.6(7) 


105.0(3) 


44.4(2) 


19.11(7) 


8.30(7) 


ALPGEN 


3271(1) 


717.4(5) 


267.4(4) 


105.4(2) 


43.7(2) 


18.68(8) 


7.88(5) 


AMEGIC 


3270(1) 


717.3(7) 


266.3(6) 


105.4(3) 


44.3(5) 







a [pb] 


Number of jets 


77 + jets 





1 


2 


3 


4 


5 


6 


Comix 


45.64(5) 


25.23(6) 


18.57(6) 


9.64(4) 


4.65(2) 


2.07(2) 


0.88(3) 


AMEGIC 


45.66(3) 


25.41(6) 


18.81(7) 


9.82(3) 









Tab. 8 Cross sections a in the MC4LHC comparison [32] setup. In parentheses the statistical error is stated 
in units of the last digit of the cross section. Note that for AMEGIC++ and COMIX all subprocesses are 
considered, while ALPGEN is restricted to up to four quarks. 
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a [nb] 


Number of jets 


7 + jets 


1 


2 


3 


4 


5 


6 


Comix 
AMEGIC 


89.5(2) 
89.6(1) 


19.65(6) 
19.60(5) 


7.52(3) 
7.59(2) 


2.664(8) 
2.64(2) 


1.000(5) 


0.387(2) 




u [pb] 


Number of jets 


e~v e + bb + jets 





1 


2 


3 


4 


5 


Comix 

ALPGEN 

AMEGIC 


9.40(2) 
9.34(4) 
9.37(1) 


9.81(3) 
9.85(6) 
9.86(2) 


6.82(5) 
6.82(6) 
6.98(3) 


4.32(4) 
4.18(7) 
4.31(6) 


2.47(2) 
2.39(5) 


1.28(2) 




a [pb] 


Number of jets 


e~e + + bb + jets 





1 


2 


3 


4 


5 


Comix 

ALPGEN 

AMEGIC 


18.90(3) 
18.95(8) 
18.90(2) 


6.81(2) 
6.80(3) 
6.82(2) 


3.07(3) 
2.97(2) 
3.06(4) 


1.536(9) 
1.501(9) 


0.763(6) 
0.78(1) 


0.37(1) 



Tab. 9 Cross sections a in the MC4LHC comparison 32; setup. In parentheses the statistical error is stated 
in units of the last digit of the cross section. Note that for AMEGIC++ and COMIX all subprocesses 
are considered, while ALPGEN is restricted to up to four quarks. 





nb] 




Number of jets n 


QCD jets 


7 


8 


55 


-> ng 




49.1(4) 


14.2(3) 


99 


— > (n — 


2)5 2q 


17.0(1) 


6.0(1) 


99 


— > (n — 


4)g4q 


1.69(1) 


0.74(5) 


99 


— > {n — 


6)569 


0.0401(5) 


0.0297(8) 


99 


-^8q 






0.000158(5) 


99 


—* (n — 


1)9 H 


30.5(2) 


9.9(2) 


99 


— > (n — 


3)5 3q 


8.46(6) 


3.38(6) 


gq 


— > (n — 


5)g5q 


0.565(7) 


0.332(8) 


91 


—t{n — 


7)9 7q 


0.00501(6) 


0.0067(2) 


qq 


— > ng 




0.0209(1) 


0.0067(1) 


qq 


— > (n — 


2)g2q 


4.97(4) 


1.84(3) 


qq 


— > (n — 


4)g 4q 


1.044(9) 


0.477(9) 


qq 


— > (n — 


6)5 6q 


0.0374(3) 


0.0291(5) 


qq 


->8 9 






0.000223(4) 



a [pb] 


Number of jets n 


e + v e 


+ QCD jets 


5 


6 


qq 






0.256(2) 


0.0768(6) 


qq -> 


e + ^ e (n — 


2)5 2-7 


6.49(3) 


2.92(3) 


qq -> 


e + v e (n — 


4)54-7 


0.591(3) 


0.449(8) 


qq -> 


e + v e 6-7 






0.00640(7) 


59 


e + v e (n — 


1)51-7 


20.0(1) 


8.21(8) 


59 


e + v e (n — 


3)5 3q 


4.03(2) 


2.14(2) 


99 -> 


e + v e (n — 


5)5 5q 


0.0741(4) 


0.094(1) 


55 


e + i> e (n — 


2)g2q 


2.13(1) 


0.775(5) 


55 


e + v e (n — 


4)5 4-7 


0.1817(9) 


0.1058(7) 


.93 


e + i> e 6q 






0.001403(7) 



Tab. 10 Subprocess cross sections a in the MC4LHC comparison [32] setup. In parentheses the statistical error 
is stated in units of the last digit of the cross section. 



efficiency 


Number of jets 


jets 


2 


3 


4 


5 


6 


7 


8 


e = 10~ 3 
e = 10~ 6 


9.3-10- 2 
3.1-10- 2 


7.8-10- 3 
3.8-10- 3 


2.M0- 3 

1.5-nr 3 


7.0-10- 4 
4.3-10~ 4 


3.6-10- 4 
2.4-10~ 4 


1.3-10~ 4 
9.9-10- 5 


6.M0~ 5 
5.8-10- 5 




efficiency 


Number of jets 


e + v e + jets 





1 


2 


3 


4 


5 


6 


e = 10- 3 
e = 10" 6 


1.5- 10- 1 

1.6- 10- 2 


2.4- 10" 2 

4.5- 10- 3 


9.i-nr 3 

3.3-10" 3 


2.0-10~ 3 
1.2-10" 3 


6.7-10- 4 
4.3-10" 4 


1.9-10~ 4 
1.3-10~ 4 


3.1-10- 5 
2.8-10~ 5 



Tab. 11 Efficiencies for processes in the MC4LHC comparison [32] setup. 



23 



References 



[1] G. Ossola, C. G. Papadopoulos and R. Pittau, Reducing full one-loop amplitudes to scalar integrals at 
the integrand level, Nucl. Phys. B763 (2007), [HOES |hep-ph/0609007| . R. K. Ellis, W. T. Giele and 
Z. Kunszt, A Numerical Unitarity Formalism for Evaluating One-Loop Amplitudes, JHEP 03 (2008), 
10031 |arXiv : 0708 . 2398 [hep-ph] | . W. T. Giele, Z. Kunszt and K. Melnikov, Full one-loop ampli- 
tudes from tree amplitudes, JHEP 04 (2008), E13 |arXiv: 0801 .2237 [hep-ph] | . G. Ossola, C. G. 
Papadopoulos and R. Pittau, On the Rational Terms of the one-loop amplitudes, JHEP 05 (2008), 004, 
arXi v: 0802. 1876 [hep-ph]| . |S. Catani, T. Gleisberg, F. Krauss, G. Rodrigo and J.-C. Winter[ From 
loops to trees by-passing Feynman's theorem, arXiv: 0804. 3170 [hep-ph], 

[2] G. Ossola, C. G. Papadopoulos and R. Pittau, CutTools: a program implementing the OPP reduction 
method to compute one-loop amplitudes, JHEP 03 (2008). [Q42| |arXiv : 0711 . 3596 [hep-ph] | . C. F. 
Berger et al., An Automated Implementation of On-Shell Methods for One- Loop Amplitudes, Phys. Rev. 
D78 (2008) . 10360031 |arXiv: 0803 .4180 [hep-ph] | . |W. T. Giele and G. Zanderighi[ On the Numerical 
Evaluation of One-Loop Amplitudes: the Gluonic Case, arXiv : 0805 . 2152 [hep-ph] 

[3] R. Britto, F. Cachazo and B. Feng, New Recursion Relations for Tree Amplitudes of Gluons, Nucl. Phys. 
B715 (2005), [HOIS |hep-th/0412308| . R. Britto, F. Cachazo, B. Feng andE. Witten, Direct proof of 
tree-level recursion relation in Yang-Mills theory, Phys. Rev. Lett. 94 (2005), 181602, |hep-th/0501052|. 

[4] S. D. Badger, E. W. N. Glover, V. V. Khoze and P. Svrcek, Recursion relations for gauge theory 
amplitudes with massive particles, JHEP 07 (2005). [Q25| |hep-th/0504159| . S. D. Badger, E. W. N. 
Glover and V. V. Khoze, Recursion relations for gauge theory amplitudes with massive vector bosons 
and fermions, JHEP 01 (2006), EMI |hep-th/0507161| . K. J. Ozeren and W. J. Stirling, Scatter- 
ing amplitudes with massive fermions using BCFW recursion, Eur. Phys. J. C48 (2006), 159-168, 
|hep-ph/0603071| . 

[5] F. Cachazo, P. Svrcek and E. Witten, MHV Vertices and Tree Amplitudes in Gauge Theory, JHEP 
09 ( 2004), EOll |hep-th/0403047 |. K. Risager, A direct proof of the CSW rules, JHEP 12 (2005), 
[0031 |hep -th/0508206j. P. Mansfield, The Lagrangian origin of MHV rules, JHEP 03 (2006), [037] 
[hep-th /0511264]. 

[6] S. D. Badger, E. W. N. Glover and V. V. Khoze, MHV rules for Higgs plus multi-parton amplitudes, 
JHEP 03 (2005), EH |hep-th/0412275| . T. G. Birthwright, E. W. N. Glover, V. V. Khoze and 
P. Marquard, Collinear limits in QCD from MHV rules, J HEP 07 (2005), IU6"8"l |hep-ph/0505219| . 
C. Duhr and F. Maltoni, Antenna functions from MHV rules, arXiv: 0808. 3319 [hep-ph] . 
[7] T. Gleisberg, S. Hoche, F. Krauss and R. Matyskiewicz , How to calculate colourful cross sections effi- 
ciently, arXiv : 0808 . 3672 [hep-ph] . 
[8] A. Kanaki and C. G. Papadopoulos, HELAC: a package to compute electroweak helicity amplitudes, 
Comput. Phys. Commun. 132 (2000). [SO^STB] |hep-ph/000"2082] . C. G. Papadopoulos, PHEGAS: A 
phase space generator for automatic cross-section computation, Comput. Phys. Commun. 137 (2001), 
1247-2541 |hep-ph/0007335| . |A. Cafarella, C. G. Papadopoulo s and M. Worek[ Helac-Phegas: a gener- 
ator for all parton level processes, |arXiv*l 0710 . 2427 [hep-ph] 
[9] F. Krauss, R. Kuhn and G. Soff, AMEGIC++ 1.0: A Matrix Element Generator In C++, JHEP 02 
(2002), Ell |hep-ph/0 109036| . 
[10] F. Maltoni and T. Stelzer, MadEvent: Automatic event generation with MadGraph, JHEP 02 (2003), 
E27t |hep-ph /0208156 ]. J. Alwall et al., Ma dGraph/MadEvent v+: The New Web Generation, JHEP 09 
(2007). I028 | |arXiv: 0706 .2334 [hep-ph]] . 
[11] M. L. Mangano, M. Moretti, F. Piccinini, R. Pittau and A. D. Polosa, ALPGEN, a generator for hard 

multiparton processes in hadronic collisions, JHEP 07 (2003), 001, hep-ph/0206293 . 
[12] F. A. Berends and W. T. Giele, Recursive calculations for processes with n gluons, Nucl. Phys. B306 

(1988) , [75S 

[13] F. A. Berends and W. Giele, The Six Gluon Process as an Example of Weyl-Van Der Waerden Spinor 
Calculus, Nucl. Phys. B294 (1987), 700, R. Kleiss and H. Kuijf, Multi-gluon cross-sections and five jet 
production at hadron colliders, Nucl. Phys. B312 (1989), 616, F. A. Berends, W. T. Giele and H. Kuijf, 
Exact Expressions for Processes Involving a Vector Boson and Up to Five Partons, Nucl. Phys. B321 

(1989) , 39, F. A. Berends, W. T. Giele and H. Kuijf, On six jet production at hadron colliders, Phys. 
Lett. B232 (1989), [2661 

[14] M. Dinsdale, M. Ternick and S. Weinzierl, A comparison of efficient methods for the computation of 
Born gluon amplitudes, JHEP 03 (2006), 056, |hep-ph/0602204 . 



24 



C. Duhr, S. Hoche and F. Maltoni, Color-dressed recursive relations for multi-parton amplitudes, JHEP 
08 (2006), EH |hep-ph/0607057| . 

P. D. Draggiotis, R. H. P. Kleiss and C. G. Papadopoulos, Multi-jet production in hadron collisions, 
Eur. Phys. J. C24 (2002), EOS |hep-ph/020220l] . 

F. Caravaglios and M. Moretti, An algorithm to compute Born scattering amplitudes without Feynman 
graphs, Phys. Lett. B358 (1995), EUOSS |hep-ph/9507237| . 

E. Byckling and K. Kajantie, N-particle phase space in terms of invariant momentum transfers, Nucl. 
Phys. B9 (1969),ESH3IS 

A. van Hameren and C. G. Papadopoulos, A hierarchical phase space generator for QCD antenna 
structures, Eur. Phys. J. C25 (2002), [563321 |hep-ph/0204055| . 

M. L. Mangano, S. J. Parke and Z. Xu, Duality and multi-gluon scattering, Nucl. Phys. B298 (1988), 
653. 

V. Del Duca, L. J. Dixon and F. Maltoni, New Color Decompositions for Gauge Amplitudes at Tree and 
Loop Level, Nucl. Phys. B571 (2000), [SOi |hep-ph/9910563| . V. del Duca, A. Frizzo and F. Maltoni, 
Factorization of tree QCD amplitudes in the high-energy limit and in the collinear limit, Nucl. Phys. 
B568 (2000), [HOH |hep-ph/9909464] . 

F. Maltoni, K. Paul, T. Stelzer and S. Willenbrock, Color-flow decomposition of QCD amplitudes, Phys. 
Rev. D67 (2003). IUT2U2"6l |hep-ph/020927l] . 

IL. J. Dixon! Calculating scattering amplitudes efficiently, hep-ph/9601359. 

K. Hagiwara and D. Zeppenfeld, Helicity Amplitudes for Heavy Lepton Production in e + e~ Annihilation, 
Nucl. Phys. B274 (1986), HJ 

S. Dittmaier, Weyl-van-der-Waerden formalism for helicity amplitudes of massive particles, Phys. Rev. 
D59 (1999), 016003 |hep-ph/9805445| . 

R. Kleiss and R. Pittau, Weight optimization in multichannel Monte Carlo, Comput. Phys. Commun. 
83 (1994), EOS |hep-ph/9405257| . 

G. P. Lepage, VEGAS - An Adaptive Multi- dimensional Integration Program, CLNS-80/447. 
F. James, Monte-Carlo phase space, CERN-68-15. 

F. Caravaglios, M. L. Mangano, M. Moretti and R. Pittau, A new approach to multi-jet calculations in 



hadron collisions, Nucl. Phys. B539 (1999), 215-232, hep-ph/9807570 



R. Kleiss, W. J. Stirling and S. D. Ellis, A new Monte Carlo treatment of multiparticle phase space at 
high energies, Comput. Phys. Commun. 40 (1986), 359 



G. C. Blazey et al. , Run II jet physics, hep-ex/0005012 



http : //mlm . home . cern. ch/mlm/mcwshop03/mcwshop .html 

S. Jadach, Foam: Multi- dimensional general purpose Monte Carlo generator with self-adapting simplical 
grid, Comput. Phys. Commun. 130 (2000). [244^2591 |physics/9910004| . 

S. Hoche, F. Krauss, S. Schumann and F. Siegert, A comprehensive approach to CKKW merging, in 
preparation. 

http : // www . scotgrid . ac . uk 



25 



